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Abstract. Hypervolume indicator is a commonly accepted quality mea- 
sure for comparing Pareto approximation set generated by multi-objective 
optimizers. The best known algorithm to calculate it for n points in d- 
dimensional space has a run time of 0{n d ^ 2 ) with special data structures. 
This paper presents a recursive, vertex-splitting algorithm for calculating 
the hypervolume indicator of a set of n non-comparable points in d > 2 
dimensions. It splits out multiple child hyper-cuboids which can not be 
dominated by a splitting reference point. In special, the splitting refer- 
ence point is carefully chosen to minimize the number of points in the 
child hyper-cuboids. The complexity analysis shows that the proposed 
algorithm achieves 0((|) n ) time and 0(dn 2 ) space complexity in the 
worst case. 



1 Introduction 

Optimization for multiple conflicting objectives results in more than one optimal 
solutions (known as Pareto-optimal solutions). Although one of these solutions 
is to be chosen at the end, the recent trend in evolutionary and classical multi- 
objective optimization studies have focused on approximating the set of Pareto- 
optimal solutions. However, to assess the quality of Pareto approximation set, 
special measures are needed [1]. 

Hypervolume indicator is a commonly accepted quality measure for compar- 
ing approximation set generated by multi-objective optimizers. The indicator 
measures the hypervolume of the dominated portion of the objective space by 
Pareto approximation set and has received more and more attention in recent 
years [2,3,1,4]. 

There have been some studies that discuss the issue of fast hypervolume cal- 
culation [5-8]. These algorithms partition the covered space into many cuboid- 
shaped regions, within which the approach considering the dominated hyper- 
volume as a special case of Klee's measure problem is regarded as the cur- 
rent best one. This approach [8] adopts orthogonal partition tree which requires 
0{n d / 2 ) storage and streaming variant [9]. Conceptual simplification of the im- 
plementation are concerned and thus the algorithm achieves an upper bound of 



<3(nlogn + n d / 2 ) for the hypervolume calculation. Ignoring the running time of 
sorting the points according to the d-th dimension, O(nlogn), the running time 
of this approach is exponential of the dimension of space d. 

This paper develops novel heuristics for the calculation of hypervolume in- 
dicator. Special technologies are applied and the novel approach yields upper 
bound of 0((f )") runtime and consumes 0(dn 2 ) storage. The paper is orga- 
nized as follows. In the next section, the hypervolume indicator is defined, and 
some background on its calculation is provided. Then, an algorithm is proposed 
which uses the so-called vertex-splitting technology to reduce the hypervolume. 
The complexities of the proposed algorithm are analyzed in Section 4. The last 
section concludes this paper with an open problem. 

2 Background 

Without loss of generality, for multi-objective optimization problems, if the d 
objective functions / = (/i, ••-,/<<) are considered with /j to be minimized, 
not one optimal solution but a set of good compromise solutions are obtained 
since that the objectives are commonly conflicting. The compromise solutions are 
commonly called Pareto approximation solutions and the set of them is called the 
Pareto approximation set. For a Pareto approximation set M — {3/1, y 2 , ■ • ■ , Dn} 
produced in a run of a multi-objective optimizer, where yi — (yn, . . . ,y i( i) E 
M C R d , all the solutions are non-comparable following the well-known concept 
of Pareto dominance. Specially, we say that j/j dominates yk at the j-th dimension 
if Vij < Vkj- 

The unary hypervolume indicator of a set M consists of the measure of the 
region which is simultaneously dominated by M and bounded above by a ref- 
erence point r — (ri, . . . , rd) G R d such that rj > maxj = i ; ... ;n {yij}. In the 
context of hypervolume indicator, we call the solutions in M as the domina- 
tive points. As illustrated in Fig. 1(a), the shading region consists of an or- 
thogonal polytope, and may be seen as the union of three axis-aligned hyper- 
rectangles with one common vertex, i.e., the reference point r. Another example 
in three dimensional space is shown in Fig. 1(b), where five dominative points, 
3/1 = (1,2, 3), 3/2 = (4, 3, 2), 3/3 = (5, 1,4), y 4 = (3, 5,1), 3/5 = (2,2,2.5), and the 
reference point r — (6, 6, 6) are considered. The volume is the union of the vol- 
umes of all the cuboids each of which is bounded by a vertex, where the common 
regions are counted only once. If a point yk is dominated by another point 3/3, 
the cuboid bounded by yk is completely covered by the cuboid bounded by yi. 
And thus only the non-dominated points contribute to the hypervolume. 

3 The proposed algorithm 

In other works, e.g. the work of Beume and Rudolph [8], the hyper-cuboid in 
(i-dimcnsional space are partitioned into child hyper-cuboids along the <i-th di- 
mension and then all these child hypervolumes are gathered together by the 
inclusion-exclusion principle [10]. 



(a) A hypervolume indicator in the two- 
objective case 



(b) A hypervolume indicator in the 
three-objective case. To lay out the 
cuboids well, the axes are rotated where 
the reference point is shaded 



In this paper, we step in another way. The hyper-cuboid is partitioned into 
child hyper-cuboids at some splitting reference points and then all the child hy- 
pervolumes are gathered directly. More detailed, given a point £ M, each of 
other points in M must dominated yi at some dimensions for the non-comparable 
relation. If the parts over yi are handled, the problem of calculating the hyper- 
volume bounded by M and the reference point is figured out. The additional 
part partitioned out at the j'-th dimension is also a d-dimensional hyper-cuboid 
whose vertices are ones beyond j/j at such dimension. Their projections on the 
hyperplanc orthogonal to dimension j are all dominated by yi, and thus are 
free from consideration. It should be noted that the reference point of child 
hyper-cuboid is altered to r' — (ri, . . . , yij, . . . , r^), namely the j-th coordinate 
is replaced by y^ . The other child hyper-cuboids are handled in the similar way. 
In these processes, the given point is called the splitting reference point. 

Obviously, the hyper-cuboids with more dominative points require more run 
time to calculate the hypervolumes. To reduce the whole run time for calculating 
all these child hyper-cuboids, the splitting reference point should be carefully 
selected. The strategy adopted in this paper is described as follows. 

(1) Let k = n — 1 and choose a point with the least dimensions on which the 
point dominated by other k points. 

(2) If some points tie, update k as k — 1 and then within these points, choose 
a point with the least dimensions on which the point dominated by other k 
points. 

(3) Repeat the similar process until only single point is left or k — 1. And if 
k = 1 and several points are left, the first found point is selected. 



By the above principle, as an example, not yi or other points but y$ is chosen 
as the first splitting reference point for the case shown in Fig. 1(b). Two child 
cuboids each bounded by one points and another child cuboid bounded by two 
points are generated by splitting along y 5 . This is the optimal strategy in such 
case. 

The algorithm to calculate the hypervolume is shown in Algorithm 1. Some 
major parameters are as follows. 

— int [n] [d] order The orders of all the dominative points at each dimension 
are represented by a two-dimensional array of integer. 

— int split The index of the point at which the hyper-cuboid is cut to 
generate multiple child hyper-cuboids is called split. 

— int[n] splitCount The numbers of k present in the split-th row of the 
array order are saved in splitCount, where k = 0, . . . , n — 1. 

— int[n] coveredCount The numbers of k present in the current checked 
row of the array order are save in coveredCount, where k = 0, . . . , n — 1. 

Moreover, some conventions are explained as follows. 

— The subscript of y^ begins with 1 while the index of array begins with 0. 
Thus yij is same as y[i — — 1]. 

— Assume a and b are two arrays and n is an element. a[] <= n means setting 
each element of a as n, while a[] <= 6[] means copying all the elements of b 
to a pairwise. 

— Assume S is a set and x is an element. £ <= S + {x} means appending a 
copy of x to S. 

The inputs of the algorithm are a set of non-dominated (dominative) points 
and a reference point, thus the hyper-cuboids are represented implicitly. 

In fact, when the hyper-cuboid is cut into two child hyper-cuboids, there may 
be some points dominated by the splitting reference point in the bigger cuboid, 
and thus such points could be removed from the points set H . In the proposed 
algorithm, it does not matter whether those points are removed or not. 

4 Complexity Analysis 

Before discussing the time-space complexity of the proposed algorithm, some 
properties are presented firstly. 

Lemma 1. Let Sij be the number of points dominating y% at the j-th dimension. 
Then 

(1) For d>2 and each i G {1, . . . ,n), J2j=i $ij >n—l. 

(2) Ford>2 and each j G {1, . . . , d}, J27=i % < f ( n - !)■ 

(3) For d~2 and each i G {1, . . . ,n}, $Zj=i % = n — 1. 

U) (»-!)■ 

(5) Ford>2 and each i G {1, . . . ,n}, Y?j=\ S ij < U n ~ X )- 



Algorithm 1 Calculate Hypervolume, CalcVolume(H) 

Input: The hyper-cuboid H defined by the dominative points {2/1, J/2, ■ ■ ■ ,y n } where 
Hi = (yn, ■ ■ ■ , Hid), and the reference point r — (n, . . . , r<j), namely H = 
{j/i, . . . ,y n ,r}. The initial number n of dominative points can be obtained from 
the length of H and the dimension d is known too. 

Output: The hypervolume of H, volume. 
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/* initialization */ 
if n = 1 then 

return ]J d . =1 \r 3 - yij \; 
end if 

volume <= 0; 
splitCount[] n; 

/* count the numbers of points dominating every point at each dimension */ 
for j = 1 to d do 

sort yij ,...,y nj ; 

for all i such that 1 < i < n do 

order [i — l][j — 1] <= number of points dominating ytj strictly; 

end for 
end for 

/* estimate split based on the statistical results of order */ 
for i = 1 to n do 

cover edCount[] <= 0; 

for j = 1 to d do 

cover ■edCount[order[i — 1, j — 1]]++; 

end for 

for k — n — 1 downto do 

if cover edC ount[k\ < splitC 'ount[k] then 
split •<= i; 

splitCount[] <= cover edCount[]; 
break; 
end if 
end for 
end for 

/* cut H at each dimension through the point indexed by split */ 
for j = 1 to d do 

if order[split — l][j — 1] > then 
H2 <= {}; 

for all yi in H\{y apUt , r} do 

if yij is dominated strictly by y sp m.j then 
H2 <= H2 + {y t }; 
yij y spin,]', 
end if 

/* Here yi can be removed from H if yt is dominated strictly by y av ut 
end for 
r2 r; 

r2 \j ~ !] Vs P ut,j\ 
H2^H2 + {r2}; 

volume ■<= volume + CalcVolume(H2); 
end if 
end for 

volume •<= volume + YYj=i \ r i ~ Vsplit,j\; 
return volume] 



Proof. It is clear that (2) => (4) => (5). The follows show (1), (2) and (3). 



(1) (By Contradiction.) Assume to the contrary there is some i G {1, . . . , n}, 
Ej=i % < n — 1. If this is the case, there are at least one yk where k ^ i 
such that each yij dominates y^j for all j G { 1 , . . . , d} . It follows that yi 
dominates yk, which contradicts our assumption that all the elements in 
{yi, . . . , y n } are non-comparable. 

(2) Given j, sort all y^ where i = l,...,n and label each y^ a sequence 
number I(i) which ranges from to n — 1. Thus J(i) = ^(n — 1). 
There are two cases to consider. Firstly, if all are different each other, 
then — It follows that E"=i $ij = §( n ~ !)• Secondly, if there 
are same elements within {yij, . ■ . , y n j}, without loss of generality, suppose 
Vij — Vkj and J(fc) = i"(i) + 1. Then Sij = Skj = < I{k), it follows that 
J27=i % < §( n ~ -*-)• This completes the proof. 

(3) (By contradiction.) For any y i} J2j=i^ij < n — 1 is excluded by (1) of 
this lemma. Thus Ej=i % > ^ — 1 for some is considered. If this is 
the case, we obtain Y^i=i Ej=i $ij > n ( n ~ 1)' contradicting (2) of this 
lemma, which implies £)" =1 Ej=i <**j = Z)j=i Z)"=i % ^ n ( n _ !)> namely 

Er=iEj=i *«<»(»-!)■ 

Lemma 2. Let Wj(fc) &e f/ie amount of k in all Sij where j = 1, . . . , d, namely 
= \{j '■ Sij = k ,j = 1,. . . Then 

(1) < Wj(fc) < d for any i and k; 

( 2 ) E"=i < rf /° r an 2/ fc ; 

(3) Efco fc ^( fc ) < i(« - 1) f° r an y *■ 

Proof. By the definition of 0Ji(k), it is clear that all statements follows Lemma 1. 

Lemma 3. Let f(n, d) be the runtime of Algorithm 1 to compute a hypervolume 
with n dominative points in a d-dimensional space. Then 

(1) f(n,d) + f(m,d) > f(n — l,d) + f(m+ l,d) where n — m > 1; 

(2) f(n, d) > f(m, d) + f(n — m, d) where n > m; 

(3) f(n, d) is minimal when Ej=i $ij =n—l and \5ij — Stk\ < 1 for any j and 

k; 

(4) f(n, d) is maximal when Ej=i % = |( n — 1) f or an U * an d w i(^) = % f or 
any i and each k = 0, . . . , n — 1. 

Proof. (1) and (2) are clear. 
(3) By the process of Algorithm 1, given some i, 

d 

f(n,d) = dnlogn + ^2f(Sij,d) (1) 

By (1) of Lemma 1, Ej=i > n — 1. It is clear that for a given i, it is 
necessary that Ej=i &j = ft — 1 to minimize /(n, d). In addition, all the Sij 



must share alike, i.e. \Sij — Sik\ < 1 for any j and k. If this is not the truth, 
suppose Sij — 5ik > 1. Thus by (1) of this lemma, 

f(5ij,d) + f(S ik , d) > f(Sij - 1, d) + f(5 ik + 1, d) (2) 

Let Siji = S.^ — 1 and 8ik> = Sik + 1. and <5,fc< can be modified in the 
similar way until \Sij — Sik\ < 1. This completes the proof. 
(4) By (5) of Lemma 1, Y^j=i — f ( n ~ -0- ^ is clear that for a given i, it is 

necessary that = |( n — 1) to maximize f(n,d). 

Hence Eqn. (1) is written as follows, 

n-l 

f(n,d) = dnlogn + J2"(k)f(k,d) (3) 
fe=i 

Suppose ?/j is the splitting reference point chosen by Algorithm 1, Ui(n— 1) < 
^, or else contradicting Y^i=i w i( n — ^) — ^- To maximize /(n, d) in Eqn. (3), 
let u>i(n — 1) = -. Similarly, we get Wj(n — 2) = -, . . ., ^i(l) = -, and so 
on. It is exactly X)fc=o = f ( n — !)• This completes the proof. 



4.1 Bounds of runtime at special cases 

First of all, it is clear that f(l,d) — d. By (3) of Lemma 3, the algorithm performs 
best when each Sij shares alike for the chosen i. If d > n — 1, 5^ < 1 for any j. 
Thus 

f(n,d)=dnlogn+(n-l)f(l,d) (4) 

which implies f(n,d) = Q(dn\ogn). If d < n — 1, Sij > 1 for any j. In the rough, 
we get 

TL — 1 ft 

f(n, d) = dnlogn + d ■ /( — - — , d) < dnlogn + d ■ f( — ,d) (5) 

d d 

It can be obtained from Eqn. (5) that f(n,d) = ( dn log nlog d n) even when 
f(^,d) is relaxed to f{%d). 

Fredman and Weide [11] have shown that Klee's measure problem has a lower 
bound of Q{n\ogn) for arbitrary d > 1. Just as Beume and Rudolph [8] have 
mentioned, although it is unknown what the lower bound for calculating the 
hypervolume is, it is definitely not harder than solving KMP because it is a 
special case of KMP. Therefore, there is a gap between the lower bound of the 
proposed algorithm and the actual lower bound of calculating the hypervolume. 

In the average cases, suppose that for the given splitting reference point yj, 
J2j=i % — f ( n ~ !)■ Meanwhile, each S^ shares alike, i.e. 5^ = Thus, 

TL — 1 Tl 

f(n,d) = dnlogn + d- /(— — ,d) < dnlogn + d ■ f(-,d) (6) 



which implies the runtime of the proposed algorithm is 0(dn logd ) at the given 
cases. 



4.2 Upper bound of runtime 

By (2) of Lemma 3, f(n - 1) > f(n - 1 - k) + f{k) for any k = 1, . . . , And 
by (4) of Lemma 3, at the worst cases, we have 

f(n, d) = dnlogn + i (f(n - 1) + /(n - 2) + . . . + /(2) + /(l)) 

<dnlogn+|(l + 2f^)/(n-l) (7) 
< dnlogn + §/(n - 1) 

which implies that the proposed algorithm for computing the hypcrvolumc bounded 
by n points and a reference point in (i-dimensional space has a runtime of 
at the worst cases. 

4.3 Space complexity 

Let g(n, d) be the used storage by Algorithm 1. In the proposed algorithm, every 
child hypervolume is calculated one by one. Since the storage can be reused 
after the former computation has been completed, g(n, d) is only related to the 
maximum usage of all the computations of child hypervolumes. Hence, 

g(n, d) — dn+ max {g(Sij,d) : < Sij < n — 1} (8) 

iE{l,...,n},je{l, ■■■,<!} 

Thus the upper bound of space is as follows. 

g(n,d) — dn + g(n — l,d) (9) 

where g(l) — d. It is easy to obtain an 0(dn 2 ) space upper bound for the 
proposed algorithm. 

Combining the above analyses together, we obtain the time-space complexity 
of the proposed algorithm. 

Theorem 1. The hypervolume of a hyper-cuboid bounded by n non- comparable 
points and a reference point in d-dimensional space can be computed in time 
0((|)") using 0(dn 2 ) storage. 

5 Conclusions 

A fast algorithm to calculate the hypervolume indicator of Pareto approxima- 
tion set is proposed. In the novel algorithm, the hyper-cuboid bounded by non- 
comparable points and the reference point is partitioned into many child hyper- 
cuboids along the carefully chosen splitting reference point at each dimension. 
The proposed approach is very different to the technique used in other works 
where the whole (i-dimensional volume is calculated by computing the (d — 1)- 
dimensional volume along the dimension d. Such difference results in very differ- 
ent time bounds, namely 0((^) n ) for our work and 0(jv* ) for the best previous 
result. Neither kind of technique can exceed the other completely and each has 
his strong point. Additionally, the amount of storage used by our algorithm is 



only 0(dn 2 ) even no special technique is developed to reduce the space complex- 
ity. 

As the context has mentioned, it is very important to choose appropriate 
splitting reference point for our algorithm. Well selected point can reduce number 
of points in separated parts and thus cut down the whole runtime. We do not 
know whether the strategy adopted in this paper is optimal or near optimal. 
Further investigations should be worked on. 
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